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ABSTRACT 

We investigate the gravitational wave (GW) signal generated by a population of double 
neutron-star binaries (DNS) with eccentric orbits caused by kicks during supernova 
collapse and binary evolution. The DNS population of a standard Milky-Way type 
galaxy has been studied as a function of star formation history, initial mass function 
(IMF) and metallicity and of the binary-star common-envelope ejection process. The 
model provides birth rates, merger rates and total numbers of DNS as a function of 
time. The GW signal produced by this population has been computed and expressed 
in terms of a hypothetical space GW detector (eLISA) by calculating the number of 
discrete GW signals at different confidence levels, where ‘signal’ refers to detectable 
GW strain in a given frequency-resolution element. In terms of the parameter space 
explored, the number of DNS-originating GW signals is greatest in regions of recent 
star formation, and is significant ly incr eased if metallicity is reduced from 0.02 to 
0.001, consistent with iBelczvnski et alJ (I2010T) . Increasing the IMF power-law index 
(from -2.5 to -1.5) increases the number of GW signals by a large factor. This number 
is also much higher for models where the common-envelope ejection is treated using 
the a—mechanism (energy conservation) than when using the 7—mechanism (angular- 
momentum conservation). We have estimated the total number of detectable DNS GW 
signals from the Galaxy by combining contributions from thin disc, thick disc, bulge 
and halo. The most probable numbers for an eLISA-type experiment are 0—1600 
signals per year at S/N^l, 0—900 signals per year at S/N>3, and 0—570 at S/N^5, 
coming from about 0—65, 0—60 and 0—50 resolved DNS respectively. 

Key words: Gravitational waves - neutron stars - Stars: binaries: close - Galaxy: 
structure - Galaxy: stellar content 


1 INTRODUCTION 

Most stars are members of binary or multiple star systems. 
Over 70% of massive stars (O-type stars) have a nearby 
companion which will affect their evolution, with over one 
half doing so before they lea ve the main sequence (MS) 
dSana et al .11201 j iLa.neeifeOl 2f) . If both members of a binary 
are massive enough to end their evolution as core-collapse 
supernovaiQ. after one or two explosions, the final prod¬ 
uct could become a double neutron-star binary (DNS), a 
neutron-star plus black-hole binary, or a double black-hole 
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1 Neutron stars can also be formed by accretion of Oxygen-Neon 
white dwarfs 


binary. Evidence that such systems do form is provided by 
the double pulsa r PSR J0737-3039A/B feurgav et al.l120031 : 
iLvne et al.ll200-'ll : iKiamer fc Stairsll2008l) . Such compact bi¬ 
naries are expected to be a significant source of gravitational 
wave (GW) radiation (iBarish Weisslll999l : iRicci &; Brilletl 
Il997l : IPostnov fc Yungelsonl 20061 ). They are amongst the 
sources most likely to be detected by a gravitational wave 
detector in the frequency range 10~ 5 — 100 Hz. Aasi et 
al. (2014) measured upper limits of the GW strain ampli¬ 
tudes from hundreds of pulsars using data from recent runs 
of the ground-based GW observatories - LIGO, Virgo and 
GEO600, and showed that there are good prospects for de¬ 
tections in the 10 - 1000 Hz range with the advanced LIGO 
and Virgo detectors. 


Double neutron-star and black-hole binaries play an 
important role in testing the theory of General Relativ- 
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Table 1 . Orbital periods (P 0 rb) and frequencies (/ or b) 5 eccentric¬ 
ities (e), GW timescales (tqw) and whether masses h ave been 
establ ished for known and suspected DNS binaries; after lLorimerl 

boos!') . 


DNS 

PSR 

-^orb 

(d) 

e 

M 

lo Sl0 
r Gwy r 1 

/orb 

Hz 

J0737-3039 

0.102 

0.09 

Yes 

7.9 

1.1 X 10“ 4 

J1906+0746 

0.17 

0.09 

Yes 

8.5 

6.8 x 10“ 5 

B1913+16 

0.3 

0.62 

Yes 

8.5 

3.9 x 10“ 5 

B2127+11C 

0.3 

0.68 

Yes 

8.3 

3.9 x 10“ 5 

J1756-2251 

0.32 

0.18 

Yes 

10.2 

3.6 x 10“ 5 

B1534+12 

0.4 

0.27 

Yes 

9.4 

2.9 x 10“ 5 

J1829+2456 

1.18 

0.14 

No 

10.8 

9.8 x 10“ 6 

.11518+4904 

8.6 

0.25 

No 

12.4 

1.3 x 10“ 6 

J1811-1736 

18.8 

0.83 

Yes 

13.0 

6.2 x 10“ 7 

B1820-11 

357.8 

0.79 

No 

15.8 

3.2 x 10“ 8 

In order of discovery: PSR B1913+16: Hulse &: Tavloil Jl975ll. 

PSR B1820- 

11: ILvne 

& McKenna! IU989fl. PSR B1534+12: 

IWolszczanl d 199lT), PSR B2127+11C: Prince et al.l ( 

19911"). PSR 

J1518+4904: iNice et al.l il996f). PSR J1811-1736: 

Lyne et alj 

(2000|), PSR J0737-3039A/B: iBurgav et al.l (i2003h: iLvne et al.l 

I2004I'). PSR J1829+2456: IciiamDion et alJ 120041'). PSR 

J1756-2251: 

Faulkner et al. 

(2005 

), and PSR J1906+0746: 


Lorimer et al.l ( 

20061'). 



double compact objects using the binary-star population- 
synthesis (BSPS) method. They concluded that only a few 
(2-4) NS-NS binaries in the Galaxy would have been de¬ 
tectable by the cancelled space observatory LISA. However, 
approximations for (i) the calculation of the GW signal from 
individual DNS binaries, and (ii) the employment of crucial 
initial conditions and the treatment of important physical 
processes in the BSPS method, may result in quite large 
uncertainties. 

In order to understand the signal detected by suffi¬ 
ciently sensitive GW detectors, it is necessary to characterize 
the radiation from all GW-emitting populations, including 
DNS. This paper presents a study of the GW signal from 
the Galactic population of steady-state DNSfl including the 
Galactic star-formation history, the initial-mass function, 
metallicity, and the physics of common-envelope evolution. 
It examines how this DNS population differs from the Galac¬ 
tic double-white-dwarf (DWD) population, and establishes 
a basis for computing the GW signal due to extragalactic 
DNS populations. We describe the methods used to model 
the DNS population in the Galactic disc, the emission and 
superposition of GW signals, and the reduction of the data 
in terms of a conceptual GW experiment (eLISA) in §[2] 
Major results are presented in §[BJ Implications, observations 
and previous work are discussed in §[4] The main conclusions 
are reviewed in §[5] 


ity, whilst double pulsars provide probes of magnetospheric 
physics. The accurate timing of pulsars orbiting a black hole 
can be used to constrain the strain amplitude of gravita- 
tional waves and th e physical properties of the black hole. 
(iKramer et al.ll2004l ~). 

At present (2014) 7 DNSs, have been confirmed and 3 
more are suspected (Table [lj. Of these, half have merger 
timescales shorter than a Hubble time. Half also show ec¬ 
centric orbits, even at relatively short periods. Bayesian sta¬ 
tistical analyses based on these observations indicate that an 
optimistic Galactic DNS merger rate may be up to 1.8 10~ 4 
yr~ 4 , implyin g that their number should be approximately 
a few million ( Kalogera et al . 200ll . l2004ll if we assume the 
age of the Galaxy is ~ 14 Gyr. This number is roughly 1-2 
orders of magnitude h igher than estimated from binary star 
population synthesis rtNelemans et al.l [200 ll ; lOslowski et al.l 
l201ll ; iDomin ik ct al .1 12012! ). although the uncertainty in the 
Bayesian analysis can exceed one order of magnitude. 

Theoretical investigations of the G W signal from Galac ¬ 
tic DNSs have be en carrie d o ut by (e .q.)|A]len__et__alJ |j999|); 
IBelczvnski et al.l (l2010al l ; iRosadol (120111 1: I Allen et all 
J2012I') . The stochastic GW background produced by DNS 
mergers in lo w-redshift galaxies (up to z^5) has been in¬ 
vestigated bj?^^ffimbau)&Ide^^eita7£achec3 (l2006l l (also 
see iRosa dc] (120111) ). who found that the signal should be 
detectable by the new gene ration of ground-based inter¬ 
ferometers. IZhu et all (120131 ) studied the GW background 
from compact binary mergers, and showed that, below 100 
Hz, the background depends primarily on the local merger 
rate and the average chirp mass and is independent of the 
chirp mass distribution. In addition, the effects of cosmic 
star formation rates and delay times between the forma¬ 
tion and m erger of binaries are lin ear below 100 Hz in 
their model. IBelczvnski et al.! (1201081 1 studied the GW back¬ 
ground and foreground signal from a Galactic population of 


2 METHODS 

2.1 Binary-star population synthesis 

In this section, we review the important physical processes 
in binary star evolution and the initial and boundary condi¬ 
tions for population synthesis to obtain a sample of double 
compact objects. Population synthesis, including individ¬ 
ual stellar-evolution tracks and initial conditions, was car¬ 
ried o ut u sing the method describe d bv lYu fe Jeffervl (l2010l . 
l201lll and lHurlev et all (I2000l . l2002ll in an initial study of the 
present Galactic double degenerates population. Note that 
we take some standard parameters for the Galactic disc (or 
the Galaxy) to represent a Milky-way type galaxy. 


2.1.1 Common-envelope evolution 

When one star in a binary fills its Roche lobe either by 
evolutionary expansion or by orbital shrinkage, Roche lobe 
overflow (RLOF) occurs. The Roche radius of the primary 
is given by 

R Ll 0.49 g ; /3 

a 0.6 ql /3 + ln(l + q{ /3 ) 

where a is more generally the semimajor axis of the orbit and 
qi(= mi/m 2 ) i s the mass rat i o of prima ry and secondary 
( EggIetonlll983l j . iHurlev et al] d2000l . l2002h have shown that 
there is a critical mass ratio q c of a binary star which can be 
used to distinguish the stable RLOF and common envelope 


2 A DNS merger would be clearly indicated by a strong GW 
signal at rapidly rising frequencies >0.1 Hz. The most optimistic 
Galactic DNS merger rates are < 10~ 3 yr —1 - see §03 as a function 
of several key parameters, 
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(CE) phases, where q c is a function of primary mass mi, its 
core mass mi c , and the mass-transfer efficiency of the donor. 
We adopt 


<?c = 


f 1.67 — x + 2 


( mic 

V mi 


/2-13, 


( 2 ) 


where x = 0.3 is the exponent of the mass-radius relation 
at co nstant luminosity for giant stars ll Hurley et al.l I 2000 L 
1200211 . 

The calculation of the orbital parameters ( e.g. orbital 
separation) of a binary after CE ejection in our model is 
based on one of two assumptions: either 1 ) angular momen¬ 
tum conservation ( 7 —mechanism) or 2 ) energy conservation 
(a—mechanism). 

In the first case, we consider the angular momentum lost 
by a binary system undergoing non-conservative mass trans¬ 
fer to be de scribed by the decrease of primary mass times 
a fac tor 7 (IPaczynski fc Ziolkowskil 119671 : iNelemans et al.l 

l200dk 


Ji - Jf 
Ji 


m 1 — mic 

T-* 

mi + m 2 


(3) 


where J) is the orbital angular momentum of the pre-mass 
transfer binary; Jf is the final orbital angular momentum 
after CE ejection; mi and mi c is the primary mass and its 
core mass respectively; m 2 is the secondary mass. 

Combining the above equation with the fraction of an¬ 
gular momentum lost during the mass transfer, Ji — Jf, and 
Kepler’s law, we have the ratio of final to initial orbital sep¬ 
aration 


Of 

Oi 


( m 1 \ 2 /mic +m 2 
\mic ) \ M 


) (l 7 


— mic \ 2 

M ) ’ 


(4) 


where oi and Of are the orbital separations before and after 
the CE phase; M = mi + m 2 is the su m of the primary 
and se condary mass before the CE phase. INelemans fe Touu 
(l2005li investigated the mass-transfer phase of the progeni¬ 
tors of white dwarfs in binaries employing the 7 -mechanism 
based on 10 observed systems and deduced a value of 7 in 
the range of 1.4 - 1.7. In order to investigate the influence 
of angular momentum loss on the rates of DCOs, we here 
adopt 7 =1.3 and 1.5. 

In the second case, CE ejection of a binary star requires 
that the envelope binding energy, including gravitational- 
binding and recombination energies, must represent a sig¬ 
nificant fraction of the orbital energy (I Webbinklll984l 'l. 


G(mi - mic)mi _ /Gmi c m 2 Gmim 2 \ 


Ar Ll 


V 2 df 


2a ; J ’ 


(5) 


where A is a structure parameter depending on the evolu¬ 
tionary state of the donor, oce is the CE ejection efficiency 
representing how much orbital energy was used to eject the 
CE, rLf is the Roche lobe radius of the donor at th e onset of 
mass transfer, and G is the gravitational constant (IWebbinkl 
1 19841 '). Rearranging in the form of Eq[f] we have 

1 2 (mi - micjai 

aAm 2 rL! 


at 

_ mi c 

d[ 

7771 



In this paper, we adopt A = 1.0, and a = 0.5 and 1.0. 

A major difference between the two CE ejection for¬ 
mulations is that energy conservation implies a signihcant 


spiral-in stage in order to eject the envelope (if we only con¬ 
sider the orbital energy as the main engine) whilst the 7 — 
mechanism does not, implying that the orbital separation 
can be larger after CE ejection in the latter case. In partic¬ 
ular, under the assumption of no external moment imposed 
on a conservative mass-transfer binary, the angular momen¬ 
tum of the binary must be conservative, and the final orbital 
separation can be written as 

«f I"(mi — Am)(m 2 + Am)l 2 

5 V * / 

a\ [ 77117712 

where Am is the fraction of mass transferred from the pri¬ 
mary to the secondary. This means that in conservative evo¬ 
lution, if mi > m 2 prior to mass transfer, the orbital sepa¬ 
ration incfreases after mass transfer. 

Althoug h both CE ejection formulations can repr oduce 
observations llNelemans fe Toutll2005l : IWebbinkll2008ll via a 
variation of free parameters, considering both conservation 
laws may be a better approach to constrain the CE evolution 
and final orbit of a binary. Note that in this paper we neglect 
viscosity, friction between the CE and the stellar cores, and 
the potential nuclear and chemical energy in the system. 


2.1.2 Gravitational radiation, magnetic braking, and tidal 
interaction 

Other mechanisms which reduce the orbital separation of a 
binary system include gravitational radiation and magnetic 
braking. A close compact binary system driven by gravi¬ 
tational radiation may eventually undergo a mass transfer 
phase, ultimately leading to coalescence. Using the average 
energy (E) and angular momentum (J or b) loss during one 
orbital period, we deduce the decay of orbital separation and 
eccentricity with respect to time to be 

da 64 G 3 mim 2 (mi + m 2 ) 1 + |§e 2 + §|e 4 
d t = ^ (1 -e 2 ) 7 / 2 ’ 

de _ G 3 mim 2 (mi + m 2 ) + ^-e 3 

d t c 5 a 4 (1 — e 2 ) 5 / 2 

This calculation is consistent with our calculation of the GW 
signal from DNS described in 5 12.21 

Gravitational radiation could explain the formation of 
cataclysmic variables (CVs) with orbital periods less than 
3h, while magnetic braking of the tidally coupled primary 
by its own magnetic wind would account for orbita l angular- 
momentum loss from CVs with periods up to 10 h (iFaulkneil 
Il97ll : IZangrilli et al.lll997l ') . We use the formula for the rate 
of angular-momentum loss due to magnetic brakin g derived 
bv lRappaport et all dl983h and lSkumanichl (Il972f) : 

j mb = -5.83 x lO" 16 ^ ( rtJspin ) M 0 R|yr- 2 , (10) 

m V R ©y r / 

where r, m env and m are the radius, envelope mass and mass 
of a star with a convective envelope, and w S pin is the spin 
angular velocity of the star. 

Tidal interaction caused by the gravity differential plays 
an important role in the synchronization of stellar rota¬ 
tion and orbital motion, and the circularization of the or¬ 
bit. Relatively co mplete descriptions of t he tidal evolution 
have been given bv lHurlev et alj d2002l l and lBelczvnski et al.l 


( 8 ) 

(9) 
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d2008l 'l. In this paper, we adopt the same formulae and 
procedures to deal with the t i dal e v olution as b y [Hut] 
dl98lli: fZ ahnl rtl 977t): ICamobelll ll 19841 ) ; iRasio et, al.1 (|l996l ) 
and Hurley et al. ( 200211 . 


2.1.3 Formation of double neutron stars 

Our simulations assume three routes for the formation of 
neutron stars: i) if a star has core mass of m c < 2.25 Mg at 
shell helium ignition, it evolve through double-shell thermal- 
pulses up the asymptotic giant branch. The star may become 
a neutron star if its core mass grows and eventually exceeds 
2.25 Mq; ii) if the core mass of the star on the thermal- 
pulsing asymptotic giant branch does not exceed 2.25 Mq 
but it is heavy enough (m c >1.6 Mq) to become an electron- 
degenerate oxygen-neon white dwarf which may become a 
neutron star via accretion-induced collapse; iii) if a star has 
a core mass of m c >2.25 Mq at the start of the early asymp¬ 
totic giant (or red giant) branch, it will become a neutron 
star without ascending the thermal-pulsing asymptotic giant 
branch. If the core mass of a star at the time of supernova 
explosion is sufficiently high (> 7Mq), it will most likely 
become a black hole unless signifi cant mass loss takes place. 
These criteria are consistent with lHurlev et ahl d200ol 'l. 

The gravitational mass of neutron stars is calculated by 

m ns = 1.17 + 0.09m Q , (11) 

where m a represents either the mass of the carbon-oxygen 
core at the time of supernova explosion or the mass of the 
oxygen-neon core, estimated by m a = maxjmch, 0.773 m c — 
0.35} with mch being the Chandrasekhar mass. Since m a is 
in the range of ~ 1.4 — 7Mq, the masses of neutron stars 
are in the range of 1.3-1.8 Mq. This i s consistent with ob¬ 
servat ional and theoretical constraints. lLattimer &: Prakasbl 
■:. 2llllTl i show that about 83% of observed neutron stars have 
mass in the range 1 — 2 Mq, while 100% of observed neu¬ 
tron stars have mass in the range 0.8 — 2.5 Mq. For this 
pap er, we assume the r a dius of a neutron star to be 10 
km dLattimer fe Prakashl d2007t) give empirical values in the 
range 9-15 km). 

For the formation of DNS, we assume seven evolution 
channels: 

I: CE ejection + CE ejection; 

II: Stable RLOF + CE ejection; 

III: CE ejection + stable RLOF; 

IV: Stable RLOF + stable RLOF; 

V: Exposed core + CE ejection; 

VI: Solo CE ejection; 

VII: Solo RLOF. 

A binary with mass ratio q = mi/ m 2 less than some 
critical value q c will experience dynamically stable mass 
transfer if the primary fills its Roche lobe while the star 
is in the Hertzsprung gap or on the red giant branch. The 
primary will become a compact object and the orbital sep¬ 
aration will change as 

— d In a = 2d In m 2 + 2an.LOFd In mi + d ln(mi + m 2 ) (12) 


where orlof is the mass- transfer efficiency for stable Roche- 
lobe overflow (R LOF) ( Han et al.| 19951). H ere, we take 
qrlo f = 0.5 (IPaczvhski fe Ziolkowskil 19671 : iRefsdal et al.l 
Il974 ). Subsequently, if the secondary fills its Roche lobe 


while it is in the Hertzsprung gap or on the red giant branch, 
then RLOF will occur. 

If the adiabatic response of the radius of the mass 
donor is less than the change of its Roche lobe radius 
with respect to a change of mass, i.e. ( a in j^° not J < 

I a ! n ~iw RLOF I , mass transfer will be unstable and a com- 
^sinM donor y RLOF ’ 

mon envelope (CE) will form. Interaction (friction) between 
the compact cores and the CE will convert orbital energy 
into kinetic energy, heating and expanding the CE. If the 
energy conversion mechanism is sufficiently efficient, the CE 
will be expelled and a compact binary with a short orbital 
period will result (see 1) 12.1.11) . 

The above formation channels for DNSs represent the 
combination of RLOF and CE processes. Note that channel 
V occurs when the envelope of a massive primary is removed 
by a stellar wind rather than a first CE ejection. CE ejection 
following evolution of the secondary may also give rise to a 
compact binary. 


2.1.4 Neutron star kicks 


Observations of proper motions indicate that pulsars have an 
extraord inary natal veloc i ty higher than their nominal pro¬ 
genit ors |Mkrkgwgkiil97d;|Lyn^ L^ne_^ - Lorimeij 

ieyl 


I 1994 I ; lHansen fe Phinnevl 1997 : Fryer fe Kaloeera 1997). 
This may result from binary evolution ( Iben fe Tutukovl 


1 19961) and asymmetric collapse and explosio n of supernovae 
( Lai et ahlll995l . l200ll : iNordhaus et aklfioi^ ). In this paper, 


we assume that both cases can contribute to the acquisition 

of kick velocities b y neutron stars. _ 

Althoug h iFr yer fe Kalogeral (Il997l ) and 
lArzoumanian et al.l 1 20021) suggested a possible bimodal 
distribution for the kick velocities, we here simply assume 
that the kick velocities have a Maxwellian dist ri butio n 
following the best estimate of lHansen fe Phinnevl (ll997h . 
taking selection effects into account, with 


d N 
Vd«k 


- (2M 1/2 ^ e ~ v i /2 ^, 


(13) 


where Uk is the kick velocity and <Tk is its dispersion, d N/N 
is the normalized number in a kick velocity bin dt>k. We 
take <7k = 190 km s -1 , so that the probable kick velocity 
Uk p ~ 268 km s -1 . As with other parameters in our popula¬ 
tion synthesis, a Monte Carlo procedure is used to generate 
the individual kick velocities for the neutron stars. Other 
parameters associated with the kick velocity (e.g. the di¬ 
rection) are assumed to follow a unif orm distribution. For 
comparison, lHansen fe Phinnevl (ll997l ) gives a kick velocity 
distribution of pulsars with a mean value of about 250 — 300 
km s 1 and a k = 190 km s -1 . Witho ut considering any se¬ 
lection effects. ILvne fe Lorimerl (Il994l) found a mean p ulsar 
kick velocity of 450±90 km s -1 . Hobbs et al.l (I2005T) sug¬ 
gested that there is lack of evidence for the bimodal distri¬ 
bution of kick velocities and a <7k = 265 km s _1 . 

A kick adds a component to the orbital velocity of neu¬ 
tron stars imposes, leading to a change of the binary orbit. 
We correct the values of orbital parameters of neutron stars 
using Kepler's laws and the binary dynamics. We have 

a( 1 — e 2 ) =| d x v 0 | 2 /( GM ), 


(14) 
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where G is gravitational constant, M the total mass of bi¬ 
nary, a the semi-major axis of the orbit, e the eccentricity, 
d the distance vector between the two stars, and v 0 the or¬ 
bital velocity of neutron star. The orbital velocity v 0 can be 
expressed as 


d 1993h and Kroupa (12001) const r ained by the observations 
of IWielen_et_ah (Il983 Popper (Il980h and the H ipparcos 
. 1l998l : l.lahreifi fc Wielenlll997l l. 


mission (Creze et al 


In this paper, we adopt the frequently used power-law 


IMF 


5 V " , GM (j-L - _L) . (15) 

The velocity of the mass centre of the binary v c relative to 
the old mass centre can be simplified to 

( M' - Ami)v c = micVk + Am' 1 ^v 0 , (16) 

where M' is the total mass of the binary before supernova 
collapse, Am© the mass loss of the (exploding) primary due 
to the supernova collapse and explosion, and m' 2 the mass 
of secondary. From Eas. ll4lfol and Kepler’s laws, the new 
orbital parameters are determined. We use the new orbital 
parameters to calculate the GW emission from DNS. Note 
that our model for pos t-kick neutron stars is con s istent with 
the dynamic model in iBrandt &: Podsiadlowskil (Il995|j and 
I Hurley et al.1 (l2002l f . 

2.1.5 Initial conditions for population synthesis 

In order to obtain a sample of compact binaries in the Galac¬ 
tic thin disc which is comparable with observations, we have 
performed a Monte-Carlo simulation in which we need the 
six physical inputs described below. In this study, we only 
investigate the effect of the first three. We use the differ¬ 
ent cases in our simulations to obtain information on the 
population of DNS and their GW radiation in other Galac¬ 
tic components via mass-scaling. The Galactic structure is 
described in § 12.1.61 and the results are presented in § 13.41 

(i) We adopt three star formation (SF) models in our bi¬ 
nary population synthesis to see the influence of SF history 
on the population of compact binaries. These are: 
Instantaneous SF: a single star burst at the formation of the 
thin disc with a constant SF rate of 132.9M© yr -1 from 0 
to 391 Myr. followed by no SF from 391 Myr to 10 Gyr; 
Constant SF: SF occurs at a constant rate of 5.2 Mq yr -1 
from 0 to 10 Gyr; 

Quasi-exponential SF: the star formation rate S is the com¬ 
bination of a major star-forming process (the first term of 
the following function) and a minor star formation (the sec¬ 
ond term of the function), 

S(t s {) = 7.92e -(tsf)/T + 0.09(t sf ) Meyr ” 1 (IT) 

where t s f is t he tim e of star formation, and r — 9 Gyr 
dYu fe Jeffervl l20icil 'l . which produces «3.5 M©yr -1 at the 
current epoch. 

All three models produce a thin-disc star mass of ~ 
5.2 x 10 10 Mq at the thin-disc age t = 10 Gyr. Obser¬ 
vations place the current thin-disk S F rate in the range 
» 3 — SMpyr - 1 dSmith et aid Il978l : iTimmes et aid 1 19971 : 
[Diehl et a.lJ i2006h and imply that Instantaneous SF is highly 
implausible for the SFH in the thin disc. It is retained here 
for comparison. 

(ii) The initial mass function (IMF) can be constrained 
by the local luminosity function, stellar density and po¬ 
tential. The IMF for th e Galactic compone n ts may be dif- 
ferent as indicated by iRobin et aid d2003h . iKroupa et aid 


£(m) = Am a , 0.1 ^ m 100M© (18) 

where m is the primary mass; ((m)dm is the number of 
stars in the mass interval mtomldm; A is the normaliza¬ 
tion coefficient determin ed by A £,(m)dm = 1. Since the 
results of lKroupa et al.l dl993ll and iKroupal d200lh indicate 
that a is in the range of —1.3 to —2.7, we take a = —1.5 
and —2.5 for comparison. 

(iii) We have adopted a metallicity Z = 0.02 (Population I) 
and 0.001 (Population II). 

(iv) We assume a constant mass-ratio distribution 
dMazeh et al .1119921 : iGoldberg fe Mazehlll994l l , 

n(l/q) = 1, 0.001 < 1/q < 1. (19) 

The inverse mass ratio has a minimum value of 0.001 since 
the maximum and minimum mass of MS stars is 100 and 
0.1 Mq in our simulations. 

(v) We employ the distri bution of in i tial or bital separations 
used bv iHanl lll998ll and lHan et al.l (l2003h . where they as¬ 
sume that all stars are members of binary systems and that 
the distribution of separations is constant in log a (a is the 
separation) for wide binaries and falls off smoothly at close 
separations: 

da _ f asep(^) ,a<ao, ^o) 

dn \ a sep , ao < a < ai. 

where a sep « 0.070, ao = 1077©, ai = 5.75 x 10®.R© = 
0.13 pc, k « 1.2. This distribution implies that the number 
of binary systems per logarithmic interval is constant. In 
addition, approximately 50 per cent of all systems are binary 
stars with orbital periods of less than 100 yr. These binaries 
are excluded when, during evolution, the condition that the 
sum of their radii exceeds their initial orbital separation is 
satisfied. 

(vi) The distr ibution of initial eccentricities of b inaries fol¬ 
lows P e — 2e (lHeggielll975l : iNelemans et al.ll200lf) . 

2.1.6 Galactic mass distribution - potential, and 
rotational velocity 

We consider the Galaxy to comprise three components, 
namely the bulge, disc, and a dark matter halo. We assume 
that the position of the Sun is given by its distance ftom the 
Galactic centre i? sun = 8.5 k pc and heig ht ab ove the Galac¬ 
tic plane z S un = 16.5 pc dFreudenreichlll998h . Our approxi¬ 
mation for the Galactic density distribution is summarized 
in Tabled The detailed expressions are described as follows. 

(1) We adopt a normal density distribution for 
the spherical bulge with a cut-off radius of 3.5 kpc 
dNelemans et al.1121)11 1:, 

Pb(r) = -^ 3 e ~ (r/ro) M©pc~ 3 , (21) 

where r is the radius from the center of the Galaxy, ro = 
0.5 kpc is the bulg e scale length, and Mb = 2.0 x 10 10 is the 
mass of bulge. IRobin et al] (120031 1 suggest that the structure 
of the inner bulge (< 1° from the Galactic center) is not 
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Table 2. Density laws and associated parameters, r is the spher¬ 
ical radius from the center of the Galaxy and ro is bulge scale 
length; R and 2 : are the natural cylindrical coordinates of the ax- 
isymmetric disc, Jir is the scale length of the disc, h z is the scale 
height of the thin disc, h' z is the scale height of the thick disc; a' 
is the radius of the halo and a' Q is a constant; p c is the central 
mass density. 



density law 

constants 

(kpc) 

Pc 

(M 0 pc -3 ) 

Bulge 

e — ( r / r o ) 2 

t-q = 0.5 

Mb = 12.73 

Thin disc 

e - R / h Rsech 2 (-z/h z ) 

h R = 2.5 

h z = 0.352 

- 1.881 

Anh^hz 

Thick disc 

e -R/ h rc e -^/ h ' z 

h R = 2.5 

h' z = 1.158 

M ,( k , = 0.0286 
4-rvhj^h^ 

Halo 

K 1 + (w -) 2 )] -1 
a 0 

o,q = 2.7 

0.108 



yet well constrained observationally. Consequently we here 
focus on the outer bulge and make no allowance for any 
additional contribution to the compact-binary population 
from the central region. 

W e use the potential proposed by I Miyamoto 8z Nagail 
(I 1975 I) in cylindrical coordinates to calculate the rotational 
velocity of stars in the bulge. 

(2) We model the thin and thick disc components of the 
Galaxy using a squared hyperbolic secant plus exponential 
distribution expressed as: 

p^ z ) = ^h e ~ R/hRp{z) M ° pc " 3 ’ ( 22 ) 

where R and z are the natural cylindrical coordinates of the 
axisymmetric disc, and Hr = 2.5 kpc is the scale length of 
the disc, h = h z for the thin disc, h = h' z for the thick disc, 
and Md = Mt n = 5.2 x 10 10 Mq is the mass of the thin disc; 
Md = M t k = 2.6 x 10 9 Mq is the mass of the thick disc. 
p(z ) is the distribution in z, with: 

p(z) = sech 2 (—z/h z ) (thin disc) (23) 

and 

p(z) = e~ z ! h% (thick disc), (24) 


where h z = 0.352 kpc is the scale height of the thin disc 
and h' z = 1.158 kpc is the scale height of the thick disc. We 
neglect the age and mass dependence of the scale height. 

The Miyamoto & Nagai potential in cylindrical coordi¬ 
nates is also used to calculate the rotational velocity of stars 
in the disk. 

(3)For the halo, we employ a rela tively simple density 
distribution which is co nsist ent with Caldwell fc Ostrikerl 
(ll98lT) , IPaczvnskil (Il990l4 and iRobin et al.l (l2003h : 


Ph(a') = p Ch 1 + Gr 


(25) 


where a' is the radius of the halo, p Ch = 0.108 Mq pc -3 
and a' 0 = 2.7 kpc. 

For the dark matter halo, we adopted the potential of 


Figure 1. Rotational velocity as a function of galactocentric 
distance R from the Galatic model, showing the contribution due 
to different components, i.e. the bulge, thin disc + thick disc, 
and halo including dark matter. Th e da s hed line indicates the 
observational estimate bv lBrand fc Blit3 ( 19931 ). The spheroidal 
component due to the interstellar medium was not considered 
separately. 


CaldweU_^_Ostrikei (1198lh . The observational estimate by 


Brand fe Blitz (1993]) is included for comparison. 


Fig. Q] demonstrates the influence of the Galactic model 
on the rotation curve of the Milky Way. 

From the Galactic model, the total mass of the halo 
including dark matter is 4.5 x 10 11 Mq inside a sphere of 
radius 50 kpc. We only focus on the baryonic mass in the halo 
which is considered to be 5 x 1 Q 10 Mq, constrained by the 
density of double white dwarfs (lYu fe Jeffenlkoich . With 
the SFR adopted here, the baryonic mass in the bulge and 
disc is at least 2 x 10 10 Mq and 5.5 x 10 10 Mq respectively, 
implying that our model assumes no dark matter component 
localized to the bulge or the thin disc. 

Combining the Galactic model and the mass of the 
Galactic components, the stellar density in the solar neigh¬ 
bourhood is 0.064Mqpc -3 , of which 6.27 x 10 -2 Mqpc -3 is 


in the thin disc, 9.4 x 10 4 Mgpc 3 is in the thick disc, and 
2.18 x 10~ 5 Mqpc~ 3 is in the halo. This is consistent with 

"~ 3 


the Flipparcos result, 0.076 ± 0.015Mqpc -3 (ICreze et al.l 
Il998i ). The local dark matter density in our model is about 
O.OlMgpc 


-3 


2.2 Gravitational waves from double neutron stars 

We calculate the GW strain amplitude from a single NS- 
NS bi nary by linearizing the e quations of general rela - 
tivity (IPeters fe Mathewa Il963l : lLandau fe Lifshit3 Il975l ). 
We assume the following properties of a binary. 1) The 
masses of the two components are mi and m 2 respectively. 
Therefore the total mass M = mi + m 2 , and the reduced 
mass p, = 2) The semi-major axis of the or¬ 

bit is a. 3) The eccentricity is e. 4) From Kepler, the or- 





















































bital separation d = i^cosCy) and the angular velocity 

p = < ' Gt ' I ^/ 2 ^p 2 ° 3 p 2 > ' > ^ , where p is the angle between the 
orbital separation and the ^-direction of an arbitrary Carte¬ 
sian coordinate in the plane of the binary orbit, and G is 
the gravitational constant. We take the z—direction perpen¬ 
dicular to the orbit plane and the origin is at the center of 
mass. Clearly, dp = / Q P ° rb pdt, where P or b is the orbital 
period. 

The c omponents of the strain amplitude can be ex¬ 
pressed as dLandau fc Lifshit zlll975r ) 


( 


^XX 

h yyi 

hzx 


hyi y 

^yy 

Hzy 


^xz 

Hyz 

hzz 


2 G 

3c 4 jRb 


■^xx -^-xy -^xz 

^yx Ayy Ay Z 

^zx A Z y Azz 


(26) 

where A a p represents the second order differential of the 
mass-quadrupole tensor with respect to time, suffix a/3 de¬ 
notes the direction, Rb is the distance form the observer, and 
c is the speed of light. Since the mass quadrupole tensor is 


( 



Ay 

Ay 

Ay 



/rd 2 (3 cos 2 p — 1) 
3fid 2 sin p cos p 


3 fid 2 cos p sin p 
fid 2 (3 sin 2 p — 1) 
0 



we have 

Ax = —C 26 (cos 2p + e cos 3 p) + Az, 

Ay = C 26 (e 2 + e cos p + cos 2p + ecos 3 p) + Az, 
Az = — C 22 (e 2 +e cosip), 

Ay = —C 26 (sin 2p + 2 e sin p cos 2 p + e sin 3 p), 
A'x = Acy, 


(28) 


C 2 


GMfi 
a( 1 -e 2 )' 


From Eas. 1261 l27l and 1281 it can be shown that the GW 
from a binary with circular orbit is monochromatic with 
frequency 2/P or b. When the orbit is eccentric, the wave be¬ 
comes polychromatic and there is a frequency broadening in 
the GW signal. 

The average power of the GW radiated from two point 
masses over one orbital period can be obtained by solv¬ 
ing the third order differential of the mass-quadrupole 
tensor with respect t o time dPeters fe Mathewsl 1 19631 : 
[Landau fe Lifshital 19751 1 . We quote the result 


Lgw = 


32 /G 4 fi 2 M c 


z(e), 


(29) 
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Tl 2 

g{n,e) = — {[J n - 2 (ne) - 2eJ„-i( ne ) + — J TO (ne) 

+ 2 eJ n+ i(ne) - J n+ 2 (ne )] 2 
+ (1 - e 2 )[J„_ 2 (ne) - 2 J n (ne) + J„ + 2 (ne )] 2 

+ 3 ^[A(ne)] 2 }, 

where J n (ne) are Bessell functions of the first kind and n = 
1,2,3,.... Since the sum of the power in each harmonic is 
equal to the total power emitted from the binary, we have 


^fg(n,e) = z(e). (33) 

n=1 

After a ma thematical transformation based on th e 
equations above dNelemans et al.ll200ll : lYu fe Jeffervll201(t ). 
we obtain the strain amplitude h(n, e ) in the vicinity of the 
Earth at GW frequency f n in the n th harmonic as 


h(n, e) = h n 

, W2(2„)=» (A 51 ) 

= 1.14 x 10 -21 


X 


g(n,e ) \ 1/2 f M\ 5/3 (P 0 r b \ 2/3 /PA 1 
n 2 ) \M q ) \ h ) \kpcj ’ 

(34) 


fn — n/Porb, (35) 

where M = fi 3 ^ 5 M 2 ^ 5 is the so-called chirp mass. Eqs. m 
13211341 and [35] are the main equations used to calculate the 
power and strain amplitude of the GW signal from one in¬ 
dividual NS-NS pair in frequency space. These equations 
also tell us that the power and strain amplitude of the GW 
signal from a binary consisting of two point-masses can be 
determined by four parameters - the chirp mass, orbital pe¬ 
riod, eccentricity and distance. In this paper, we refer to the 
first three as orbital parameters, and use the mass of each 
component instead of the chirp mass. 

The energy flux of GW waves can be expressed as 


P = 


c 3 » 2 u 2 

16t rG 



(36) 


where fl = 2nf is the angular frequency. We define the spec¬ 
tral function of the energy flux as 


1 dP _ 2 tt 2 d (f 2 h 2 ) 
~ ^df ~ 3 H 2 df 


where p c = 3Hp/8nG is the critical mass density of the 
present universe with Hp » 73 km s " 1 Mp c -1 being the 
Hubble constant ([Freedman fe Madord [20101 1 . We show the 
GW energy flux spectrum for selected cases to assist com¬ 
parison with other work which uses this metric. 


1 + (73/24)e 2 + (37/96)e 4 

(1 _ g2)7/2 


(30) 


After Fourier analysi s of Kepler motion, we obt ain the power 
in the n th harmonic dPeters fe Mathewslll963l l 


r n — 
Vgw — 


32 

5 


/G 4 fi 2 M 3 \ 


g(n,e), 


(31) 


2.3 The sensitivity of eLISA 

Evolved-LISA (eLISA) is designed as an updated version of 
LISA, a space-based GW detector, consisting of 1 mother 
and 2 daughter satellites flying in formation to form a 
Michelson-type Laser interferometer with an arm length 
of lxlO 6 km (see http://www.elisa-ngo.org/). Noise arises 
mainly from the displacement noise (including the noise 
































caused by laser tracking system and other factors) and par¬ 
asitic forces on the proof mass o f an accelerometer (acceler¬ 
ation noise) (jLarson et alJl200oH ^l We can convert the noise 
signal to an equivalent GW signal in frequency space by 

hi = (38) 

where S n is the total strain noise spectral density, h{ is the 
root sp ectral density and R is the GW transfer function 
given bv lLarson et al.l (j2000l ). In the simulations, we take the 
displacement noise to be 1.1 10 ” 11 mHz -1 / 2 at 10 mHz, and 
the acceleration noise to be 3 10 ~ 15 ms _ 2 Hz _1//2 at 10 mHz. 
By comparison, the arm length of LISA would have been 
5 10 6 km. The displacement noise and acceleration noise 
would have been 4 10 -11 mHz -1 / 2 and 3 10 ~ 15 ms _ 2 Hz -1 ^ 2 . 

For a continuous monochromatic source, such as a NS- 
NS binary with a circular orbit, which is observed over a 
time T, the root spectral density will appear in a Fourier 
spect rum as a single spectral line in the form (lLarson et al.l 

l 2000 h 

hi = hVr. ( 39 ) 

So, for an observation time T = 1 yr, the root strain ampli¬ 
tude spectral density hi = 5.62 x 10 3 4 h Hz -1 / 2 . 

To demonstrate the detectability of the predicted GW 
signal due to the galactic DNS population, we show the ex¬ 
pected eLISA sensitivity for S/N=l in the figures in §[3j 

2.4 Data reduction 

We reduce the simulated GW signal by using its mean in¬ 
tegrated value (h). To do this, we first choose a frequency 
interval A f' which is greater than the interval used to com¬ 
pute the simulations (i.e. A/ = lyr -1 ). We then calculate 
the mean value ( h) of the strain amplitude and its standard 
deviation a^ in this large frequency interval using 

Eti hi _2 ZLi(hi-(h)) 2 

(h) — -:- , a (h) — -:-> (40) 

3 3 

where j represents the number of small-frequency inter¬ 
vals in the large frequency interval. In this paper, we take 
A log f' = 0.03, so j is also a function of frequency. We plot 
(h) as a function of GW frequency in all figures unless spec¬ 
ified otherwise. In each panel, we also show the maximum 
standard deviation which represents the maximum uncer¬ 
tainty in each large frequency interval in different cases. 

2.5 Computation procedure 

In order to compute the superposition of the GW signal from 
the entire DNS population in the Galactic disc and compare 
the signal with the sensitivity of the proposed detectors, we 
need to know their birth rates, merger rates, present number, 
space, mass, eccentricity and orbital distributions, we have 
adopted the following procedure to obtain the above physical 
properties. For the Galactic thin disc with age t^ BC , having 
a total mass M tn (tdisc) = / 0 tdisc S(t a f)df s f 0: 

3 The values of parameters to calculate the total noise can be 
found on http://www.elisa-ngo.org/ 

4 We here neglect the inter stel la r medi u m wh ich makes up about 
20% of the thin disc mass iRobin et al.|[2003l l. 


Table 3. Main parameters and their values in our simulation. See 
§ 12. II for an explanation of the parameters. 





CEE(aA =) 

1.0 0.5 

CEA( 7 =) 
1.3 1.5 


SFH 

IMF(ff =) 






Con 

-1.5 

Cl 

C3 

C25 

C27 



-2.5 

C2 

C4 

C26 

C28 

Z=0.02 

Exp 

-1.5 

C5 

C7 

C29 

C31 



-2.5 

C6 

C8 

C30 

C32 


Inst 

-1.5 

C9 

Cll 

C33 

C35 



-2.5 

CIO 

C12 

C34 

C36 


SFH 

IMF(cr =) 






Con 

-1.5 

C13 

C15 

C37 

C39 



-2.5 

C14 

C16 

C38 

C40 

Z=0.001 

Exp 

-1.5 

C17 

C19 

C41 

C43 



-2.5 

C18 

C20 

C42 

C44 


Inst 

-1.5 

C21 

C23 

C45 

C47 



-2.5 

C22 

C24 

C46 

C48 


(i) For each star formation epoch in the disc, we calculate 
a sample distribution of k coeval MS binaries having a total 
mass m p and generated by the four Monte-Carlo simulation 
parameters m, q, a and e. 

(ii) We follow the evolution of each primordial binary in a 
time grid consisting of many time intervals to establish the 
properties of DNSs formed from the above MS binaries up to 
fdisc- From the timescales from MS binary formation to DNS 
formation and from DNS formation to DNS merger, we can 
obtain the contribution function from each star formation 
epoch to the number of new-born and merged DNSs in a 
given time interval, as well as the total number of alive DNSs 
during that interval. 

(iii) By summing the contributions from all star forma¬ 
tion epochs in the thin disc, we can obtain all the physi¬ 
cal information of the DNS population, including birth and 
merger rates, present number, and distributions of orbital 
parameters. 

(iv) Since most of the DNSs have eccentric orbits, we 
compute the Fourier transform of each orbit to obtain the 
value of GW strain emitted at each harmonic frequency (as 
indicated by Eqs. 1321 and 1341 1. and then sort the DNSs by 
the harmonics of orbital frequency (equivalent to the GW 
frequency). 

(v) We then calculate the total strain amplitude h 2 from 
the number and distance of DNSs in each frequency bin (for 
one year observation of e-LISA, see § 12.31) . This yields a raw 
data of strain amplitude against the GW frequency. 

(vi) Finally, we use the method described in § 12.41 to re¬ 
duce the raw data, producing a reduced relation between 
the strain amplitude and GW frequency. 

In our population synthesis simulation, we started with 
16 10' primordial MS binaries, distributed equally between 
all 48 cases, and assume that tdi ac = 10 Gyr. The parameters 
in the present study are SFH (Instantaneous, Continuous 
and Quasi-exponential), IMF (a = —1.5 and —2.5), metal- 
licity (Z = 0.02 and 0.001), and two different CE evolution 
processes (a and 7 ). For each CE formalism we adopt two 
parameters (aA = 1.0 and 0.5, 7 = 1.5 and 1.3). The 48 
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Figure 2. Birth rates of DNSs in the Galactic disc in different cases. 










































































10 



0.01 0.1 1 100.01 0.1 1 10 
time (Gyr) time (Gyr) 



time (Gyr) time (Gyr) 


Figure 3. Merger rates of DNSs in the Galactic disc in different cases. 
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Table 4. Ranges in the present birth rate, merger rate and total number of DNS from binary-star population syntheses for a thin-disc 
age of 10 Gyr. Parentheses refer to simulation labels defined in Table 02 and used in Figs. 1-8. Merger rates and total numbers may span 
a significant range, even when the present birth rate is zero, because of the delay time from new-born DNS to coalescence. 


Star Formation 

CE Ejection 

Present Birth Rate 
xl0 -5 yr -1 

Present Merger Rate 
xl0 _5 yr _1 

Total Number 

Continuous 

a 

0.49 - 5.8 (C8 - C15) 

0.49 - 5.7 (C8 - C15) 

61 - 5.5 x 10 4 (C20 - C13) 

Continuous 

7 

0.31 - 27 (C32 - C37) 

0.27 - 11 (C32 - C37) 

7.1 x 10 3 - 1.9 X 10 6 (C32 - C37) 

Instantaneous 

a 

0 

0 

0 - 5.2 x 10 4 (010,12,24 - C21) 

Instantaneous 

7 

0 

0 - 2.8 (C34,36,48 - C45) 

6.4 x 10 3 - 1.6 x 10 6 (C36 - C45) 


individual cases are labelled C1-C48 and are summarized in 
TableEl 

3 RESULTS 

3.1 The double neutron-star population 

3.1.1 Birth rates, merger rates, and numbers 

The total number of DNS is a crucial input parameter for 
calculating the integrated GW signal from the Galactic DNS 
population. The evolution of the birth and merger rates of 
DNS in the thin disc are shown in Figs.[5]and Figs. [3] There 
is a significant difference in total DNS numbers between the 
instantaneous SF model and the two other SF models. Since 
we assume that the response of descendent stars to a normal¬ 
ized star-formation epoch with the same input parameters 
is an intrinsic property, there should be a quasi-linear rela¬ 
tion between the SF rate and the DNS birthrate. This can 
explain the different DNS birthrates in the three SF models. 
For example, the maximum DNS birthrate (1.46 10 _ 3 yr _1 ) 
in case C9 (instantaneous SF: SFR= 133MQyr _1 ) is about 
26 times that (5.72 lCU 5 yr -1 ) in the case Cl (constant SF: 
SFR= 5Mgyr _1 ). The IMF affects the DNS birthrates in 
a similar way. Since the total number of initial MS binaries 
in the sample is constant ( 10 7 ), a small value of the power- 
law index a (i.e. —1.5) increases the number of massive MS 
binaries and hence gives higher DNS birthrates. 

Although metallicity and CE ejection coefficients have 
no direct impact on the initial number of MS binaries, they 
can significantly affect the orbital separation of a binary 
following a CE phase. If, after first CE ejection, the orbital 
separation is sufficiently small, the binary will not survive to 
become a DNS; this consequently reduces the DNS birthrate. 
However, if a binary can survive both first and second CE 
ejection, there is an increased probability of forming a very 
short-period compact binary. From Eqs.0] and [SJ either in¬ 
creasing 7 or decreasing a can, in principle, lead to small 
orbital separation after CE ejection, and lower the DNS 
birthrate. However, in the a—formalism, this is not true for 
low metallicity (Z = 0 . 001 ) , because the envelope of stars 
with low Z have smaller radii and usually have a smaller 
energy absorption, resulting in smaller mass loss via weaker 
stellar wind and in the formation of a massive core. From 
Eq.[ 6 ] a larger core mass helps a binary to survive, and thus 
increases the probability of DNS formation. In summary, 
under the a—formalism, a lower metallicity leads to larger 
core-mass and larger orbital separation, and increases the 


DNS formation probability. Decreasing the CE ejection co¬ 
efficient a results in a smaller orbital separation and lowers 
the DNS formation probability. In the 7 —formalism, the sit¬ 
uation is not so obvious, but the orbital separation of a bi¬ 
nary after a CE phase is also associated with the core-mass 
of primary. The DNS birthrate depends on a competition 
between the metallicity and the CE ejection coefficients. 

Once a DNS has formed, its orbital separation and evo¬ 
lution is controlled by gravitational radiation and the mag¬ 
netic field. The time scale for a MS binary to become a 
DNS is generally ~ 10 — 40 Myr. The merger time de¬ 
pends on the initial properties of the DNS (i.e. mass, or¬ 
bital period, and eccentricity). Assuming the DNS orbital- 
period number distribution to be constant in logarithm, 
(dN/da oc a -1 ), and the lifetime of an individual DNS to 
be t oc a 4 (i.e. da/dt oc a -3 ), then the lifetime number dis¬ 
tribution dN/dt oc t _1 . 

However, we argue that the DNS orbital-period num¬ 
ber distribution is not simply a constant in logarithm. We 
compute the number evolution of the DNS population (see 
Fig.[3]| using two methods to ensure consistency. First, we 
count the number from each star formation epoch. Second 
we calculate the integral of the difference between birth rate 
and merger rate. Our results show that the number of DNS 
in general depends on the total evolutionary mass of the 
thin disc. Small differences in the predicted number at the 
present epoch are found by using different metallicity and 
CE ejection coefficients. 

Ranges in the current birth rates, merger rates and 
numbers of Galactic DNS assuming a thin-disc age of 10 Gyr 
as given by our binary-star population-synthesis calculations 
are shown in Tabled 

We draw attention to case C12 (Z = 0.02, a\ = 0.5) in 
which, in our calculation and as a consequence of having or¬ 
bital periods < 6000 s, DNSs can be born and merge within 
the same computational time interval. This phenomenon 
leads to the total number of DNS always being zero in Fig. [5] 
(there is no blue dashed line in the panel showing C12). The 
similar cases CIO and C24 can also produce DNSs with short 
orbital periods (< 10 6 s, but longer than in C12). The DNSs 
in these two cases also merge quickly, resulting in zero DNSs 
at the present age of the thin disc (10 Gyr). For compari¬ 
son, as many as ~ 10 6 DNSs may be present at one time 
assuming a 7 —mechanism, that is 1 — 10 6 times higher than 
assuming an a—mechanism. Even in the extreme case of in¬ 
stantaneous SF, the present number can be ~ 10 4 — 10 6 . 
This indicates that the 7 —mechanism can produce DNSs 
with orbital periods much longer than a—mechanism. 
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Figure 4. The history of number of DNSs in the Galactic disc in different cases. Note that the alive number of DNSs in cases C4, C8. 
and C12 is zero so we can not see lines for them in this figure. 
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Figure 5. The distribution of chirp mass of new-born DNSs in the Galactic disc in different cases. The number distributions (y-scale) 
are normalized to the unity. 
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Figure 6. The distribution of eccentricity of new-born DNSs in the Galactic disc in different cases. The number distributions (y-scale) 
are normalized to the unity. 
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3.1.2 Chirpmass and eccentricity 

The distributions of chirp mass and eccentricity of the new¬ 
born DNS population are plotted in Figs.[5]and[6] Since our 
calculation is based on a number of coeval MS stars and 
SF has an overall effect on the calculation, the distribu¬ 
tions of the physical variables of the DNS population ( i.e. 
mass, eccentricity, and orbital period) are here only influ¬ 
enced by the metallicity, IMF and CE ejection coefficient. 
Since the masses of neutron stars described in Eq.[TT]in our 
calculation are in the range 1.3 — 1.8 Mg, we expect their 
chirp masses (see Eg. 1341) to lie in the range 1.1 — 1.6 Mg, 
as shown. Our results show that, in general, a lower Z or 
higher IMF index increases the fraction of DNS with higher 
chirp mass. Lower Z leads to less mass loss and eventually 
leaves a more massive core, and while a higher IMF index 
leads to more massive MS stars. In addition, we find that 
lower a or higher 7 tend to give a more concentrated distri¬ 
bution at smaller chirp mass. With no tidal interactions, the 
eccentricity distribution of the new-born DNS population is 
centered around 0.6 — 0.8 assuming an a—mechanism and 
around 0.5 — 0.7 assuming a 7 —mechanism. However, the 
eccentric orbits of the DNSs are likely to be circularized by 
gravitational wave radiation and/or tidal interaction during 
the common-envelope phase and subsequent evolution. 

3.1.3 Orbital periods 

The difference in the orbital period distributions due to 
the two CE ejection models is shown in Fig. [7] For the 
a—mechanism, most new-born DNSs have short orbital pe¬ 
riods with a peak in the range 200 — 10 000 s. Increasing the 
IMF index and/or Z gives more DNSs with longer orbital 
periods. The shortest orbital period in this model is about 
200s (a = 0.5, cr = —2.5), while the longest can be more 
than 10 8 s (a = 1.0, a = —1.5, Z = 0.001). As argued al¬ 
ready, the 7 —mechanism yields DNSs with longer orbital 
periods. The shortest orbital period in this model is about 
3 000 s. In many cases the new-born DNSs can have orbital 
periods longer than 10 8 s. These results are consistent with 
the values of birth and merger rates in Figs. [2] and E] 

3.1.4 Formation channels 

The fractions of DNSs from different binary-star formation 
channnels are summarized for each model in Table [7] For a- 
mechanism models it is found that most short-period DNSs 
(Fo r b < 10 4 s) come from the CE+CE channel, while longer 
period DNSs generally come from channels with at least 
once stable RLOF process. The dramatic decrease in the 
number of DNSs with short orbital periods (< 10 4 s) in 
the 7 -mechanism models is most likely due to the decrease 
of number of double CE ejection events. We find that the 
number fraction of DNSs with short orbital periods may be 
up to 10% in the two 7 -mechanism models where double-CE 
events contribute over 80% of the DNSs. 

3.1.5 The influence of kicks of neutron stars 

As initialized, the kick velocities of neutron stars in our sim¬ 
ulations obey a Maxwellian distribution with a maximum 
likelihood (in 3D) being ~ 268 km s -1 . With a standard 


Table 5. The fraction of DNSs from different formation channels, 
expressed as a percentage, assuming constant star-formation. 


CE 

Ejection 

Model 

z 

<T 

CE+CE 

% 

RLOF+CE 

CE+RLOF 

% 

other 

% 

Cl 

a 

1.0 

0.02 

- 1.5 

87.2 

1.6 

11.2 

C2 

1.0 

0.02 

- 2.5 

82.1 

0.0 

17.9 

C3 

0.5 

0.02 

- 1.5 

93.6 

1.0 

5.4 

C4 

0.5 

0.02 

- 2.5 

100.0 

0.0 

0.0 

C13 

1.0 

0.001 

- 1.5 

58.0 

15.9 

26.1 

C14 

1.0 

0.001 

- 2.5 

58.1 

21.0 

21.0 

C15 

0.5 

0.001 

- 1.5 

95.3 

1.6 

3.1 

C16 

0.5 

0.001 

- 2.5 

100.0 

0.0 

0.0 

C25 

7 

1.3 

0.02 

- 1.5 

44.4 

55.5 

0.1 

C26 

1.3 

0.02 

- 2.5 

48.1 

51.9 

0.0 

C27 

1.5 

0.02 

- 1.5 

94.3 

4.6 

1.1 

C28 

1.5 

0.02 

- 2.5 

85.7 

0.0 

14.3 

C37 

1.3 

0.001 

- 1.5 

37.9 

60.4 

1.7 

C38 

1.3 

0.001 

- 2.5 

47.7 

52.3 

0.0 

C39 

1.5 

0.001 

- 1.5 

64.7 

33.8 

2.5 

C40 

1.5 

0.001 

- 2.5 

59.5 

35.1 

5.4 


deviation of 190 km s - , kick velocities can lie in the range 
10 - 1200 km s -1 . Since the new orbit of a binary is derived 
from the kick velocity IS I2.1.4II . it has a crucial influence on 
the distribution of eccentricity and orbital period. Kick ve¬ 
locities mean that about 95% of DNSs have eccentric orbits, 
with about 90% from supernova kicks, and about 5% from 
binary evolution. Our results indicate that young DNSs may 
have highly eccentric orbits, although gravitational wave ra¬ 
diation (and magnetic braking) tends to circularize the or¬ 
bits during late evolution. By c omparison, all the b inary 
white dwarfs in the simulations of lYu fe Jeffervl J201(t ) have 
e = 0 due to the tidal interaction and zero kick velocity. 

Considering the extreme case of the neutron star ve¬ 
locity to be the sum of the maximum kick velocity and the 
Galactic rotation velocity, i.e. 1200+250 km s _1 , their max¬ 
imum speed is about 1.5 kpc Myr - ; the average will be 
closer to 268+250 km s -1 = 0.54 kpc Myr -1 . The measured 
GW strain depends on the distance; the change in strain over 
the course of one year due to neutron star kick velocities will 
be negligible. 

3.2 The GW strain-amplitude — frequency 
relation 

Fig. [ 8 ] shows the DNS GW signal for each model in Tablc[3] 
as a strain- amplitude frequency relation, reduced using 
the method described in §E31 so as to simplify the informa¬ 
tion contained. The probability to detect a GW signal from 
a DNS population using e-LISA, assuming an instantaneous 
SF model, is very low or close to to zero. If the SF rate 
is continuous or there is significant star formation in the 
recent past, we expect a high probability to detect a GW 
signal from a DNS population, assuming common-envelope 
ejection follows the a—mechanism. If CE ejection follows a 
7 —mechanism, DNS GW detectability with e-LISA is also 
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Figure 7. The distribution of orbital periods of new-born DNSs in the Galactic disc in different cases. The number distributions (y-scale) 
are normalized to the unity. 
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Figure 8. The gravitational wave signal from a DNS popula tion in the Galactic disc for different cases. The black line shows the 
gravitational wave noise due to double white dwarfs taken from lYu Sz Jeffervl d2010T) . There is no gravitational wave signal from double 
neutron stars (DNS) in cases C4, C8, CIO, C12 and C24 since the current total number of DNS in these cases is zero. The cyan lines 
show the sensitivity of eLISA at S/N = 1. The upper-left lines denotes the maximum deviation of the GW strain from the averaged 
value. We take the bin size of A log(//Hz)=0.03. 
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very low. The 7 -mechanism can produce more DNSs at fre¬ 
quencies, but only at frequencies about 2 orders of magni¬ 
tude lower than the best working frequency of e-LISA (10 ~ 2 
Hz). 

Most of the strain-amplitude relations show oscillations 
in ( h(f )), sometimes with amplitudes of more than 1 dex. 
The DNS chirp masses in our calculation lie in the range 
1.1 — 1.6 Mq, so can affect ( h ) by a factor of 1.87. ~95% 
of DNSs in our model lie at distances between 3.92 and 
13.92 kpc, affecting the signal by a factor of up to 3.55. These 
two parameters are unlikely to explain the oscillation. The 
distribution of DNS orbital periods should be the most im¬ 
portant factor. This can be understood as follows. 

We first assume that any probability distribution func¬ 
tion can be fitted by a normal or multi-normal distribution, 
so that the simplest distribution function is the normal dis¬ 
tribution. We apply this distribution to the orbital periods 
to understand the ( h(f )) relation. The frequency distribu¬ 
tion can be expressed as 


d N _ dN dP orb _ y? 2 d N 
df ~ dP orb ' d f - 2^ nf ' dP orb ■ 


If the distribution of P orb is normal, then 


dJV 

d7 


71 max 

nf~ 2 • C n exp 

n =1 


(nf 1 


— Por b ,C >) 2 


a 


2 

p 


(41) 


(42) 


To find the extrema / n of the function to be summed in 
Eg. 1421 we set its first derivative with respect to frequency 
/ equal to zero. After rearranging, 


f2 , nP orb 0 : - 
In H 5 In 



(43) 


where f n = n 


2(7 p 


and r/ = P OT b,o/ap. If op « 7 , 


the frequency distribution may have n max maximum val¬ 
ues, and the interval between two successive maxima is 

— »)+%/17 2 +4 


2 (7 p 


. If ap rj, the interval becomes negligible. 


The total strain amplitude h 2 in one frequency bin has 
been calculated as the sum of the strain amplitude h 2 of 
those binaries in the frequency bin, expressed as 


dN a f 

-rr A/ 

hl= Y. h '- ( 44 ) 

JV=1 


If we assume the binaries in a frequency bin have similar 
chirp mass, the above equation becomes 

hi = C h ^rf i/3 n~ 2 g(n,e), (45) 

where Ch is a constant. Combining with Eg 1421 we obtain 
an expression for the strain amplitude generated from DNS 
binaries with normally distributed orbital periods: 


,2 „ — 1 j,—2/3 ( (n/^-Porb, 0 f\ , , 

h t =C h y,n / exp I- 2 - \g(n,e). 

n=1 k a P ' 

(46) 

To find the extrema, we take the derivative d/d/ of the 
function in the sum and obtain 


r2 , „ nPorbfi n 

J nh 1 <J 9 J nh 
<T Z 



= 0 , 


(47) 


which has the solutions / n h = n 


-3r ; + A /97) 2 + 12 
2 (7 p 


. The prop¬ 


erties of / n h are similar to those of f n . In addition, when 

V '^ > 1, /hm /m* 

Consequently, if the orbital periods satisfy a perfect nor¬ 
mal distribution, an oscillation can appear in the ( h(f )) re¬ 
lation under the condition that neighboring frequencies coin¬ 
cide with extrema of (h). However, since the orbital periods 
are most likely the sum of normal distributions with different 
Por b ,o an d crp, the frequencies where the strain amplitudes 
have extrema relate to 7 and ap. 

This analysis demonstrates the following. If a popula¬ 
tion of Galactic DNS has a relatively large distribution of 
orbital period (large crp), it will generate a smooth GW 
background. However, with small ap, or a tight distribu¬ 
tion of period, we see a set of harmonics in the frequency 
range 10 -6 — 10 -2 Hz. All the DNSs in the sample are effec¬ 
tively “in tune”. With larger ap this effect is smeared out. 
This means that under some conditions ( e.g. P OT b/a 1), 
the oscillation of space caused by the GW radiation from 
a DNS population at low frequencies can be regularly en¬ 
hanced. We note that the oscillation in our simulations is 
hardly detectable by eLISA. 

For comparison, Fig. [ 8 ] also shows the reduced GW sig¬ 
nal from double white dwarf s (DWDs) simulated by the 
method in I Yu fe Jeffervl (l201(t) . The signals from DWDs oc¬ 
cupy (l$ 20 %) of the observation frequency intervals for one 
year of observation, and they should have a different polar¬ 
ization pattern from the signal from DNSs. This means that 
future observations should be able to distinguish the DNS 
and DWD signals. 

Since several authors plot gravitational energy flux 
rather than strain amplitude, Fig. [9] shows the same infor¬ 
mation as Fig. [5] for DNS in terms of the spectral energy 
distrubution ( S ), calculated from Eq.|33and reduced using 

the method in §[231 


3.3 Discrete Gravitational-Wave Signals 

In previous work (I Yu fe Jeffervl l20ld . 1201 il l we have used 
the term ‘resolved’ GW systems to refer to those binaries 
in which all of the GW strain measured within a single fre¬ 
quency interval over a given integration period arises from 
a single system. Since all of the DNSs in our simulation are 
in eccentric orbits, their GW emission may be spread over 
many frequency intervals. Thus the concept of a ‘resolved’ 
DNS may not be useful. We prefer to use the term ‘discrete 
GW signal’ to represent a frequency interval that contains 
a GW signal greater than some noise threshold over a given 
observation period. The number of such signals serves as a 
proxy to estimate what might be expected from future ob¬ 
servations. 

We hence define the number Ah CTi 3 CTj 5 CT of discrete GW 
signals exceeding the noise threshold, where the subscript 
defines the confidence level. Thus ler means that Ni a fre¬ 
quency intervals show a GW signal with S/NJST, Ns a gives 
the number of frequency intervals with S/N^3, and so on. 
The results are listed in Table[ 6 ] We find that a lower metal- 
licity and a larger IMF power-law index give a higher prob¬ 
ability for the detection of a GW signal. The highest prob¬ 
ability occurs in case C15 (Z = 0.001, a = 0.5, a = —1.5 
and constant SF rate). This case gives 57195 discrete sig- 
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Figure 9. The spectra distribution of gravitational wave (GW) flux calculated by Eq. l37l from DNS population in the Galactic disc 
in different cases. Note that there is no GW signals from double neutron stars (DNS) in cases C4, C8, CIO, C12 and C24, and the 
GW signals in cases Cll, C24, C36, C38, C42 and C46 are out of the scale of the coordinates of this figure. We take the bin size of 
Alog(//Hz) =0.03. 
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nals at la and 19933 discrete signals at 5a. The detection- 
probability for a discrete GW signal is much lower in a 
7 —mechanism model than in an equivalent a— mechanism 
model. For 7 -mechanism models, the highest probability oc¬ 
curs for C31 ( Z = 0.02, 7 = 1.5, a = —1.5 and exponential 
SF rate), with 75 signals at 1 <t and none at 5er. 

Some DNSs with large eccentricities may generate a de¬ 
tectable GW signal at multiple harmonics of the orbital fre¬ 
quency. For example, in simulation C3, a DNS with m = 
1.33+ 1.31M 0 , P OI b = 198.4 s, e = 0.508 and R h = 2.81 kpc 
generates GW signals at the fundamental frequency 0.00504 
Hz and its harmonics n = 2 — 13 {i.e. 0.01008 Hz, 0.01512 
Hz, 0.02016 Hz, ...) at S/N > 5. The numerical error in the 
frequencies is less than 0.02%. The detection of GW signals 
in harmonic series would help to constrain the DNS orbital 
parameters. Note that in addition to the signal caused by 
the DNS in this example, there is another weak GW sig¬ 
nal at frequencies of 0.02016, 0.02520, 0.04032, 0.0504 and 
0.06552 Hz (n =4, 5, 8 , 10 and 13) respectively caused by 
other DNSs. 

In order to investigate the potential number of ‘resolv¬ 
able’ DNSs, TablcQshows the numbers (N) of DNSs which 
have at least one detectable GW harmonic. This number is 
not equivalent to the number which contribute to the GW 
signals in Table[ 6 ] where a GW signal may exceed the de¬ 
tection criterion due to the combined strain of two or more 
DNS at the same frequency. Tablc[7]does not include DNS 
which generate GW signals below the detector threshold. 
Hence, for example, cases C15 and C39 show 57195 and 56 
GW signals (at lcr), representing the combined contributions 
from about 500 - 950 and 1-13 DNSs in each case. In case 
C39, no individual DNS generates a detectable GW signal. 
Figurc fTTl illustrates the number of DNSs per frequency bin 
as a function of GW frequency in cases C15 and C39. We 
finally emphasize that the GW signals in Table [ 6 ] come from 
a much smaller number of individual sources as indicated by 
the numbers of resolved sources shown in Table [7] 


3.4 GW signal from the Galactic components 

Amongst our experiments, cases C2, C4, C 6 , C 8 , C26, C28, 
C30 and C32 better represent the bulge and thin disc for 
SF history, IMF and metallicity. Howeve r, due to the lower 
total mass of the bulge ~ 1 — 2 x 10 10 tfm jRobin et af]|2003l : 
iBelczvnski et alll2010bl : lYu fc .Teffervll201ol '). the number of 
GW signals from the DNSs in the bulge should be less than 
that in the thin disc by a factor of 2.5 — 5. Since the thick 
disc has a much smaller total stellar mass (2.6 x 10 9 Mq 
dRobin et al.ll2003l : lYu fc .TeffervH201(l H. the GW signal from 
DNSs in the thick disk should be less than that in the thin 
disc by a factor of about 20. The halo has a quite different 
SF history (SF only at early epochs), IMF (a = —1.5), and 
metallicity (Z = 0 . 001 ) from the thin disc, but would have a 
similar stellar mass (5 x 10 10 Mq). Numerical experiments 
C21, C23, C45 and C47 better reflect conditions in the halo. 
Our results show that very rare GW signals from the DNSs 
in the halo can be obse rved, which support the results of 
IBelczvnski et all J2010bh . Note that the influence of distance 

may be negligible._ 

In the model of I Yu fo Jeffervl J20ld h the mean distances 
and standard deviations of stars in the bulge, the disc, and 
the halo are, respectively, 8.50 ± 0.28 kpc, 8.92 ± 2.47 kpc, 


Table 6. The number of gravitational wave signals in each sim¬ 
ulation. 
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and 8.62 ± 1.13 kpc. In fact, a realistic GW signal will lie 
somewhere between the cases studied. For example, if the 
mean metallicity of the thin disc lies between 0.001 and 0.02 
dPanter et ahllioosl ; IBelczvnski et al.ll2010bl '). a more realis¬ 
tic case would be a combination of C2 and C14, or C4 and 
C16 and so on. We here approximate the GW contribution 
from the Galaxy as a whole by correcting the thin disc con¬ 
tribution ( N ) for the bulge (+1V/2.5), thick disk (+7V/20), 

and halo (= 0 )._ _ 

IBelczvnski et all d2010bl ') concluded that the number 
of resolved DNSs in their population synthesis model 
should be abou t 6 , w hich is consistent with the results of 
iNelemans et al.l d200ll 'l . For an eLISA-type observatory, and 
assuming case C2 (Table 4) is approximately representative 
of the Galaxy as a whole, one year of observation will yield 
about 1633 observable DNS-induced GW signals, 906 signals 
at S/N=3, and 577 at S/N=5. 


4 DISCUSSION 

We have used an evolutionary model to simulate the birth 
and merger rates and the total number of DNSs in the 
Galactic thin disc. Our results on rates and total number 
are roug hl y consisten t with the observational constraints 
l|Kalogera et al.ll2004f l and with bina ry-star formation and 
evolution theory calc ulated by others llNelemans et al 1200 ll : 
iDominik et al.ll2012f) . However. INelemans et al~l ~( 2001 1 only 
give an optimistic val ue for the rates and n umber assum¬ 
ing the 7 —mechanism. iKalogera et al.l (120041 ') report results 
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Figure 10. Density map representing the gravitational wave (GW) strain amplitude of individual double neutron stars as a function 
of GW frequency in case C3. The red symbols represent the harmonic GW strain from a DNS (m = 1.33 + 1.31Mq, P OI h = 198.4 s, 
e = 0.508, and ft), = 2.81 kpc) with S/N> 5 (i.e. at 5cr level, see section [373} . The colored lines show the sensitivity of eLISA, upper line: 
S/N= 5, lower line: S/N= 1. We adopt a bin size Alog(//Hz) = 0.01 and A log h = 0.05. The grey-scale (right) shows the numbers of 
GW signals per bin. 


based on observation and statistical theory; they neglect the 
dependency of the rates and nu mber on initial condit ions 
and stellar evolution parameters. iDominik et al.l (12012 ) in¬ 
vestigated the influence of metallicity and the a—mechanism 
CE ejection coefficient, they did not investigate different SF 
and IMF models. In this paper we have systematically in¬ 
vestigated the influence of all these initial conditions and 
stellar evolution parameters on the rates and total number 
of DNSs. 

Using the binary-star population synthesis model, we 
have calculated the stable GW emission from the long-lived 
DNS population which may be observable by (for example) 
eLISA in the frequency range 0.001 — 0.1 Hz. In fact, these 
calculations can be used to estimate the GW signal from 
other Galactic components, i. e. the bulge, the thick disc, and 
the halo. The SF rates in the bulge and thin disc are most 
likely to be continuous, and both share a common metallic¬ 
ity (Z = 0.02). The power-law indices of the IMF in these 
regions are also similar, being —2.35 in the b ulge and —2.5 
in the thin disc for stars with mass m > 1 Mm (IZ occali e t al.l 
l200d : iKroupa et al JlToOOl : iKroupalboO ll : iRobin et al.ll200 J b 

iRosadol (1201 ll ) adopts an analytic approach to study 
the GW background from binary systems, and concludes 
that no background signal generated by DNSs and DWDs 


is detectable in the frequency band of ground-based GW 
detectors. Our calculations support this conclusion. 

Since the GW signal from DNSs is significantly affected 
by SF history, recent SF regions in the bulge and spiral arms 
of the Galaxy may contain many GW sources. By compar¬ 
ison, although the Large and Small Magellanic Clouds may 
have experienced bursts of star formation peaking r oughly 2, 
0.5, 0.1 and 0.012 Gvr ago_jHiarris_^_Zmitsky 20091 ) and 2.5, 
0.4, and 0.06 Gyr ago dHarris fe Zaritskvll2004 ). respectively, 
the low Magellanic-Cloud SF rate of roughly 0.2 M 0 yr -1 
implies the GW signal from Magellanic-Cloud DNSs should 
be negligible compared with that in the Galaxy. 

A caveat affecting all of our simulations is that, in prac¬ 
tice, the number of detectable DNS present in the Galaxy 
at any one time is small. Even in cases where the num¬ 
ber of GW signals exceeds 10 3 , Table [7] shows that these 
arise from a relatively small number of ‘resolvable’ DNS, 
with eccentric systems producing signals at multiple har¬ 
monics. The DNS eccentricity distribution in our model is 
essentially that defined by the common-envelope interaction 
physics, followed by tidal evolution as discussed in §2.1.2. 
While observations show that highly eccentric DNS do exist 
at long period (cf. Table 1), tidal evolution to short-period 
systems is not tested. Detailed models for the influence of 
kick velocity distribution and galactic structure on the or- 
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Figure 11. Number of double neutron stars (DNSs) per frequency bin as a function of gravitational wave frequency in cases C15 and 
C39. Note that there may be overlap of DNSs in different frequency bins. There are 56 data points in the left panel and 57195 data 
points in the right panel, corresponding to 56 and 57195 GW signals at I a level in Table [6] respectively. We adopt a bin size A/ = yr —1 
Hz. 


bital parameters of DNSs and therefore their GW signal are 
also needed to extend our model to different types of galaxy. 
More seriously, since the predicted GW signal is based on 
small numbers amongst which the eccentricity distribution 
is determined by a Monte-Carlo distribution applied to the 
parent population, uncertainties on numbers in Table 5 are 
more likely to scale as y/N/N, the numbers of DNS in Table 
6 , than as \/~N/N. 


5 CONCLUSION 

Using the methods of binary-star population synthesis, we 
have investigated the Galactic double neutron-star (DNS) 
population as a function of initial mass function, star- 
formation history, metallicity and common-envelope physics. 
The parameter space explored is larger than covered pre¬ 
viously, and includes a volume corresponding to best esti¬ 
mates of the present-day Galaxy. We have computed the 
gravitational-wave (GW) signal that would be generated by 
these theoretical populations, including the multi-frequency 
signal generated by DNSs in elliptical orbits. We have 
also explored the probability of likelihood of detecting GW 
from DNS from a conceptual space-borne GW observatory 
(eLISA). The main conclusions are: 

(1) Observable GW from double neutron stars are 
more likely in low metallicity environments, in environments 


where there is a high proportion of massive stars, and in re¬ 
gions of recent star formation. The first two statements are 
relatively obvious; low metallicity produces less luminous 
and, hence, smaller giants, so that binaries are more tightly 
bound at the point of common-envelope ejection, whilst a 
higher IMF power-law index produces more high mass stars, 
and therefore more DNS. If the galactic metallicity is re¬ 
duced from Z = 0.02 to 0.001, the number of observable 
GW signals from DNS increases by a factor up to rs 30 in 
the two continuous star formation models (e. q. C 6 and C18), 
and consistent with iBelczvnski et all (l2010al ) . The influence 
of the IMF is even stronger than that of metallicity; increas¬ 
ing the power-law index a from —2.5 to —1.5 increases the 
number of observable GW signals from DNS from 0 to 2434- 
57195 1 -a detections in the frequency range 0.0001 — 0.1 Hz. 
The fraction of DNS arising from star formation within the 
last Gyr is almost 100% of the total DNS population. 

(2) Observable GW from DNS reflect the physics of 
the common-envelope ejection mechanism. If conservation 
of energy dominates (a—formalism), DNS GW are more 
likely to be observed than if conservation of angular mo¬ 
mentum ( 7 —formalism) dominates the physics. The peak 
frequency at which the average strain amplitude has a max¬ 
imum value in these two mechanisms is also different. The 
peak frequency for the a—formalism is in the range of 0.001 
- 0.01 Hz, while the peak frequency for the 7 —formalism is 
in the range of 10 -5 —10 -3 Hz. Observation of the DNS GW 
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Table 7. The number of DNSs which have at least one detectable 
GW signal in each simulation. 


Case 

Nla 

N 3 a 

y 5ct 

Case 

N la 

y 3ct 

Nsa 

a 

Cl 

1273 

927 

743 

C13 

2434 

1446 

1031 

C2 

45 

42 

35 

C14 

1471 

812 

546 

C3 

543 

487 

429 

C15 

1473 

1291 

1115 

C4 

0 

0 

0 

C16 

38 

33 

28 

C5 

1220 

896 

729 

C17 

2005 

1194 

835 

C6 

34 

32 

26 

C18 

1049 

585 

379 

C7 

538 

474 

409 

C19 

1412 

1226 

1078 

C8 

0 

0 

0 

C20 

35 

30 

26 

C9 

0 

0 

0 

C21 

0 

0 

0 

CIO 

0 

0 

0 

C22 

0 

0 

0 

Cll 

0 

0 

0 

C23 

0 

0 

0 

C12 

0 

0 

0 

C24 

0 

0 

0 

7 

C25 

0 

0 

0 

C37 

2 

0 

0 

C26 

0 

0 

0 

C38 

1 

0 

0 

C27 

7 

1 

0 

C39 

3 

0 

0 

C28 

0 

0 

0 

C40 

0 

0 

0 

C29 

0 

0 

0 

C41 

2 

1 

0 

C30 

0 

0 

0 

C42 

2 

0 

0 

C31 

11 

1 

0 

C43 

1 

0 

0 

C32 

0 

0 

0 

C44 

0 

0 

0 

C33 

0 

0 

0 

C45 

0 

0 

0 

C34 

0 

0 

0 

C46 

0 

0 

0 

C35 

0 

0 

0 

C47 

0 

0 

0 

C36 

0 

0 

0 

C48 

0 

0 

0 


spectrum will therefore help to constrain common-envelope 
ejection physics. 

(3) Young DNSs most likely have eccentric orbits re¬ 
sulting mostly from kick velocities imparted during super¬ 
nova collapse and partly binary evolution. This creates a 
harmonic structure in the GW radiation of DNSs. 

(4) Current observations indicate the most realistic val¬ 
ues for the physical parameters in the Galactic disc corre¬ 
spond with our cases C2 and C6, i.e. Z = 0.02, a\ = 1.0, a = 
—2.5 and either constant (C2) or exponentially decreasing 
(C6) star formation. We therefore expect that one year of 
observation with a GW observatory such as eLISA will de¬ 
tect approximately 0—1600 observable GW signals caused 
by DNSs at S/N^l, 0—900 signals at S/N^3, and 0—570 
signals at S/N^5 in the Galaxy between 10~ 5 and 1 Hz, 
coming from about 0—65, 0—60 and 0—50 resolved DNS. 
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